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Abstract — Three  extensions  to  radio-frequency  (RF)  tomogra¬ 
phy  for  imaging  of  voids  under  wide  areas  of  regard  are  pre¬ 
sented.  These  extensions  are  motivated  by  three  challenges.  One 
challenge  is  the  lateral  wave,  which  propagates  in  proximity  of 
the  air-earth  interface  and  represents  the  predominant  radiation 
mechanism  for  wide-area  surveillance,  sensing  of  denied  terrain, 
or  close-in  sensing.  A  second  challenge  is  the  direct-path  coupling 
between  transmitters  (Txs)  and  receivers  (Rxs),  that  affects  the 
measurements.  A  third  challenge  is  the  generation  of  clutter  by 
the  unknown  distribution  of  anomalies  embedded  in  the  ground. 
These  challenges  are  addressed  and  solved  using  the  following 
strategies:  1)  A  forward  model  for  RF  tomography  that  accounts 
for  lateral  waves  expressed  in  closed  form  (for  fast  computation); 
2)  a  strategy  that  reduces  the  direct-path  coupling  between  any 
Tx-Rx  pair;  and  3)  an  improved  inversion  scheme  that  is  ro¬ 
bust  with  respect  to  uoise,  clutter,  and  high  attenuation.  A  finite- 
difference  time  domain  simulation  of  a  scenario  representing 
close-in  sensing  of  a  denied  area  is  performed,  and  reconstructed 
images  obtained  using  the  improved  and  the  classical  models  of  RF 
tomography  are  compared. 

Index  Terms — Green’s  functions,  ground-penetrating  radar, 
lateral  waves,  radio-frequency  tomography,  tunnel  detection. 

I.  Introduction 

The  problem  of  underground  void  detection  is  para¬ 
mount  to  secure  borders  and  sensitive  areas  and  for  search 
and  rescue  missions.  To  date,  no  underground  imaging  tech¬ 
nique  emerged  as  a  standard  for  close-in  sensing  of  wide  denied 
areas,  where  minimal  human  intervention  is  required  [1]. 

A  promising  strategy  is  introduced  in  [1]  and  [2]  where 
a  set  of  transmitters  (Txs)  and  a  set  of  receivers  (Rxs)  are 
placed  on  (or  in)  the  ground  at  arbitrary  positions.  The  Txs 
radiate  a  monochromatic  signal,  which  impinges  upon  a  buried 
dielectric  or  conductive  anomaly,  thus  generating  a  scattered 
field.  Multiple  Rxs  collect  samples  of  the  scattered  electric 
field  and  relay  this  information  to  a  base  station.  Images  of 
the  below-ground  scene  are  then  reconstructed  using  the  prin¬ 
ciples  of  RF  tomography.  The  advantages  and  mathematical 
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derivations  of  RE  tomography  for  underground  imaging  are 
discussed  in  [I].  The  approach  is  technically  valid  for  any 
sensor  disposition  and  terrain  shape  provided  that  the  Green’s 
function  characterizing  the  problem  is  properly  selected.  In  [1], 
the  Green’s  function  for  a  homogeneous  space  was  applied  due 
to  its  simplicity  of  implementation.  This  choice  has  been  proven 
to  work  satisfactorily  when  the  sensors  and  targets  are  located 
nearly  vertically  above  the  targets,  thus  avoiding  artifacts  due 
to  the  discontinuity  at  the  air-earth  interface. 

However,  practical  applications  require  wide  areas  of  investi¬ 
gations  (e.g.,  underground  networks  and  facilities),  denied  areas 
(e.g.,  sensing  of  urban  environment),  or  close-in  sensing.  In 
these  cases,  sensors  remotely  probe  underground  regions  at 
long  ranges,  and  the  propagation  of  waves  occurs  primarily 
along  the  air/ground  interface;  hence,  the  predominant  prop¬ 
agation  mode  is  the  lateral  wave  [8]-[10].  Therefore,  one 
contribution  of  this  letter  is  the  introduction  of  a  more  accurate 
forward  model  by  invoking  a  closed-form  Green’s  function  that 
accounts  for  the  air-earth  discontinuity  (see  Appendix). 

In  addition,  RE  tomography  is  based  upon  the  knowledge  of 
the  scattered  field  from  targets.  In  real  cases,  the  Rxs  are  irra¬ 
diated  by  a  strong  electromagnetic  field  due  to  the  direct  cou¬ 
pling  between  each  Tx  and  Rx  pair  (i.e.,  direct-path  coupling 
[1]).  Hence,  as  a  second  contribution,  in  Section  III,  an  efficient 
technique  that  mitigates  the  direct-path  coupling  (by  joint  Tx 
and  Rx  null  steering)  is  presented. 

Moreover,  distributed  anomalies  (e.g.,  weathered  soils)  also 
generate  a  bias  to  the  measured  scattered  field,  which  may  be 
considered  as  clutter.  A  third  contribution,  given  in  Section  IV, 
is  an  improvement  upon  the  inversion  schemes  already  dis¬ 
cussed  in  literature  [3]-[5],  based  on  the  findings  described  by 
Zhdanov  [7]  from  the  geophysical  community;  this  improved 
method  is  more  robust  with  respect  to  perturbations  (e.g., 
clutter)  of  the  measured  scattered  field. 

The  combination  of  these  three  new  strategies  improves  the 
image-reconstruction  process,  particularly  for  large  areas  of 
interest,  shallow  targets,  and  close-in  sensing,  as  shown  in 
Section  V. 

II.  Eorward  Model 

The  3-D  geometry  shown  in  Eig.  1  is  considered.  Under  the 
monochromatic  assumption  (single  work  frequency  /),  the  air 
half-space  is  modeled  as  a  free-space  medium,  while  the  ground 
half-space  is  modeled  as  a  homogeneous  medium  with  relative 
dielectric  permittivity  eo,  conductivity  au,  and  magnetic  per¬ 
meability  p,f).  The  targets  (i.e.,  tunnels,  caches,  or  voids)  are 
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Fig.  1.  Three-dimensional  geometry  for  the  model. 

assumed  to  reside  in  the  investigation  domain  D.  The  sources 
are  N  electrically  small  dipoles  (of  length  AZ)  fed  with 
current  I.  For  each  transmitting  sensor,  the  total  field  E  is 
collected  by  M  Rxs.  Both  Txs  and  Rxs  reside  inside  the  ground 
but  outside  the  investigation  domain  D.  The  relative  dielectric 
permittivity  er-(rO  ^^d  the  conductivity  cr(r')  inside  D  are  the 
unknowns  of  this  problem.  The  contrast  function  is  defined 
as[l] 


£5(r')  =  e^(r')  -60+ j  (o-(r')  -  an)  /^irfeo  (1) 

being  £0  the  free-space  dielectric  permittivity.  However,  other 
definitions  may  be  used  [3]-[5],  [8].  By  pointing  out  the  first 
order  Born  approximation  [3]-[5],  [7],  [8],  the  field  received 
by  a  dipole  oriented  along  the  direction  a^,  positioned  at 
due  to  a  transmitting  dipole  oriented  along  the  direction  a^, 
positioned  at  r^,  can  be  written  as  [1] 

E  {viyjj  =Qtxl,-G  (C,  r^)  •  +  if  (r^,  C)  +  T 

+  <3^0  JJJ  [aL  •  G(C,r')]  •  [G(r',r^) -a^]  es{r')  dr' 

D 

(2) 

where  Q  =  jcop^oAlI  for  an  electrically  small  dipole  [3];  the 
quantities  H  (multiple  scattering)  and  T  (random  noise)  rep¬ 
resent  unpredictable  perturbations  to  the  total  field.  G  is  the 
Green’s  dyadic,  to  be  chosen  according  to  the  adopted  formula¬ 
tion.  The  first  term  in  (2),  i.e.. 


Qa™  ■  G  (C,  r^)  •  a^  (3) 

describes  the  direct-path  coupling  between  a  particular  Tx  and 
Rx  pair.  The  cancellation  of  this  coupling  from  the  measured 
field  is  a  critical  problem  in  RF  tomography:  Although  it  can 
be  analytically  predicted  (and  cancelled)  using  (3),  in  practical 
cases,  its  magnitude  may  be  up  to  50-60  dB  higher  than  the 
scattered  signal.  In  these  conditions,  the  dynamic  range  of  the 
Rxs’  amplifiers  may  not  be  large  enough,  or  the  quantization 
steps  may  not  be  as  fine  as  required,  to  sample  both  scattered 
and  direct-path  coupling  fields  accordingly.  Clearly,  the  best  so¬ 
lution  is  to  cancel  the  direct-path  contribution  before  it  reaches 
the  Rx. 


Rx  dipoles  toward  the  desired  directions.  Rotation  of  dipoles 
may  be  performed  using  mechanical  devices  or  by  properly 
feeding  a  set  of  colocated  orthogonal  dipoles  (see  [1]).  The 
proposed  strategy  properly  steers  the  Tx  dipole  in  order  to 
minimize  the  field  at  the  Rx  side,  and  then  turns  the  Rx  dipole 
in  order  to  be  orthogonal  to  the  expected  direct-path  electric 
field.  Mathematically,  these  rotations  are  computed  by  solving 
a  series  of  constrained  minimization  problems  for  each  Tx  and 
Rx  pair 

minimize  ||G  (r;:;,,r^)-a*,||2  subject  to  =  1.  (4) 


In  Lagrangian  form,  it  becomes 


(aL  A)  =  (a^)  •  G^  (C,  r^)  •  G  (r^,,  r^)  •  a* 

li  -1 


-A 


m  5  ■' 
T 


(5) 


where  G^  denotes  the  Hermitian  of  G.  By  imposing 
=  0  we  obtain 

(C,r^)  •  G(r;;„r^) -a^  =  Aa^.  (6) 

Therefore,  the  a^  direction  that  minimizes  the  power  at  a 
desired  location  is  the  eigenvector  associated  with  the  smallest 
eigenvalue  of  the  matrix  G^G.  Similarly,  this  minimization 
can  be  applied  at  the  Rx  side.  Defining  the  vector 

Emi„=G(rL,r^) -a^  (7) 

as  the  electric  field  obtained  when  a^  is  chosen  according  to 
(4),  a  second  minimization  problem  can  be  formulated  as 

subject  to  (a;’„)'^  •  a^j,  =  1.  (8) 

The  minimization  is  achieved  when  is  chosen  to  be  the 
eigenvector  corresponding  to  the  smallest  eigenvalue  of  the 
outer  product  (matrix)  Emin  ■  (Emin)^-  If  the  dipole  can  be 
steered  only  over  a  horizontal  plane  (e.g.,  by  using  two  crossed 
dipoles),  the  steering  directions  are  easily  obtained  by  setting  to 
zero  the  z-component  of  each  vector  a„. 

As  tested  via  numerical  analysis,  the  application  of  these 
strategies  guarantees  an  acceptable  minimization  of  the  re¬ 
ceived  signal  due  to  direct-path  coupling.  Therefore,  the  to¬ 
tal  electric  field  can  be  reasonably  approximated  only  with 
the  first-order  scattered  contribution  from  targets  inside  the 
region  D,  i.e.,  [3] 

E(rLrm)  =L(£5(r')) 

=  III  [a;n-G(C,r')] 

D 

■  [G(r',i'U  (9) 


minimize 


•  E„ 


IV.  Inversion 


III.  Direct-Path-Coupling  Mitigation 

In  this  section,  a  direct-path  coupling  mitigation  technique  is 
introduced.  The  key  feature  is  the  possibility  to  steer  the  Tx  and 


Assuming  that  the  clutter  contribution  has  been  completely 
mitigated  and  the  field  received  by  the  sensors  is  given  by  the 
forward  model  in  (9),  the  sampled  field  at  each  Tx  and  Rx 
pair  can  be  collected  in  a  vector  e  =  {E(r^,rm)},  and  the 
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investigation  region  D  can  be  discretized  in  K  voxels,  each  one 
located  at  position  r'^,;  The  contrast  function  can  be  represented 
in  a  column  vector  e  =  {e5(r'j,)}.  After  this  discretization, 
(9)  can  be  reformulated  in  matrix  form,  i.e., 

e  =  L£  (10) 

where  L  is  generally  an  ill-conditioned  matrix. 

Several  methods  have  been  proposed  to  solve  (10),  such  as 
back  propagation  [1],  truncated  singular-value  decomposition 
(SVD)  [4],  [11],  and  Tikhonov  regularization  [7],  [11].  In  this 
letter,  we  introduce  a  refined  version  of  Tikhonov  regularization 
that  equalizes  the  sensitivity  of  each  Tx  and  Rx  pair  by  selecting 
a  proper  weighting  factor  and  introduces  a  term  that  is  able  to 
exploit  the  a  priori  information  on  the  values  of  the  dielectric 
anomalies.  Accordingly,  the  contrast  function  (as  a  function  of 
the  regularization  parameter  /3)  can  be  estimated  [7] 


Fig.  2.  Geometry  for  the  simulation.  Txs  and  Rxs  are  indicated  with  stars 
and  diamonds,  respectively,  and  the  two  tunnels  are  located  at  the  center  of  the 
scene. 


£(/?)  =  (L“W|;L  +  /3W2)  ^  ■  e  +  l3Wle°)  (11) 

where  e°  represents  the  known  dielectric  anomalies  embedded 
in  region  D,  and 

W,  =  diag(L^L)i/2  ^  diag(LL“)i/2  (12) 

are  (diagonal)  weighting  matrices  opportunely  defined  in  order 
to  minimize  the  sensitivity  of  the  system  [7]. 

In  most  cases,  the  weighting  matrices  have  a  small  dynamic 
range.  Therefore,  we  can  approximate  Wg  =  al  in  (11).  If  we 
perform  the  SVD  [11]  of  a  properly  defined  weighed  matrix 
=  W  pL  =  USV^.  (11)  becomes 

e  =  (lHL^+/3W2)-'(lH  •  e+pwl  ■  e°) 

^(VS»SV^-fa^/3I)~VVS^U^W  ^  •  e+pWl  ■  e°) 

=  Vdiag  .  e-p/JV^W^  •  e°) 

(13) 

where  Si  represents  the  dh  singular  value  of  L.^,.  The  advantage 
obtained  in  applying  (13)  is  that  the  contrast  function  as  func¬ 
tion  of  (3  can  be  computed  via  (fast)  matrix  multiplications. 

V.  Simulations  and  Conclusion 

A  simulation  is  presented  in  order  to  demonstrate  the  im¬ 
provements  achieved  by  using  the  following:  1)  lateral  waves  in 
the  forward  model;  2)  direct-path  mitigation;  and  3)  weighted 
Tikhonov  regularization.  The  test  scene  represents  a  situation 
where  sensors  are  surrounding  the  wide  (denied)  area  of  interest 
(i.e.,  close-in  sensing)  and  probe  the  region  D  mostly  via  lateral 
waves  (see  Fig.  2  for  details). 

The  targets  are  two  hollow  cylinders  (radius  1  m)  emulating 
two  tunnels,  located  with  their  axes  parallel  to  the  surface  and  at 
a  depth  z'  =  —5  m  (with  respect  to  their  center),  having  ed  =  9 
and  ao  =  5  x  10^^  S/m.  No  a  priori  information  about  dielec¬ 
tric  anomalies  in  the  scene  is  considered,  thus  =  0.  The  work 
frequency  is  5  MHz.  Placed  along  a  circle  encompassing  the 
two  tunnels  are  12  Txs  and  20  Rxs,  as  shown  in  Fig.  2.  Each 
sensor  is  emplaced  at  depth  d  =  0.25  m  beneath  the  surface. 
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Fig.  3.  Reconstructed  image  using  homogeneous  Green’s  function.  Depth 
slice:  5  m. 

The  received  electric  field  has  been  synthesized  using  the  finite- 
difference  time  domain  simulator  GPRMAX  [6]  for  each  Tx 
and  Rx  pair. 

In  the  first  simulation,  the  process  was  as  follows:  1)  The 
homogeneous  Green’s  function  [8],  [9],  having  the  properties 
of  the  soil,  is  inserted  in  the  forward  model;  2)  the  truncated 
SVD  method  [1],  [4]  is  used  for  the  inversion;  and  3)  the 
direct-path  coupling  is  assumed  to  be  completely  cancelled. 
The  reconstruction  result  is  shown  in  Fig.  3. 

In  the  second  simulation,  the  procedure  was  as  follows: 
1)  The  closed-form  Green’s  function  that  accounts  for  the  lat¬ 
eral  wave  (see  Appendix)  was  used;  2)  the  weighted  Tikhonov 
method  proposed  in  Section  IV  was  implemented;  and 
3)  the  direct-path-mitigation  algorithm  (Section  III)  was  used 
to  choose  the  direction  of  each  sensor.  The  reconstruction  result 
is  shown  in  Fig.  4. 

As  expected,  the  shallow  targets  are  mainly  irradiated  by 
the  lateral  wave  excited  at  the  air-earth  interface,  and  the 
classical  Green’s  function  for  the  homogeneous  space  cannot 
accurately  predict  the  field  value.  Conversely,  the  half-space 
Green’s  function  generates  highly  resolved  images,  particularly 
for  shallow  targets  since  it  accounts  for  the  effects  of  lateral 
wave  (see  Fig.  4).  Therefore,  we  conclude  that  when  RF  tomog¬ 
raphy  is  applied  in  wide  areas,  denied  terrain,  or  when  targets 
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Fig.  4.  Reconstructed  image  using  half-space  Green’s  function  described  in 
the  Appendix.  Depth  slice:  5  m. 

are  in  shallow  regions,  the  discussed  improvements  (i.e.,  the 
inclusion  of  the  lateral  waves  propagation  within  the  realm  of 
RF  tomography  and  its  fast  numerical  computation,  the  direct- 
path  coupling  suppression  using  steerable  radiators,  and  the 
improved  inversion  scheme  for  better  handling  the  clutter)  are 
a  suitable  choice  to  obtain  high-quality  reconstructions. 


Appendix 

Half  Space  Green’s  Function 


Half-space  Green’s  functions  have  been  proposed  in  the 
literature  [3],  [5],  [9],  [11],  but  they  are  generally  expressed 
as  asymptotic  expansions  or  in  spectral  form.  Nevertheless, 
King  et  al.  [10]  derived  explicit  closed-form  expressions  for 
the  electric  field  generated  by  horizontal  and  vertical  dipoles 
buried  in  a  lossy  medium  [under  assumption  reported  in  (16)]. 
From  King’s  formulas,  a  closed-form  expression  for  the  half¬ 
space  Green’s  function  is  derived  (valid  only  when  antennas 
and  targets  are  embedded  in  the  soil). 

A  dyadic  Green’s  function  can  be  expressed  as  follows: 
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In  this  case,  v  =  x'k  +  yy  +  zi  represents  the  observation 
point,  and  r'  =  a;'x  +  y'y  +  z'i,  represents  the  current  (phys¬ 
ical  or  equivalent)  source  position.  The  coefficients  in  (14)  are 
given  as 
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Using  the  following  notation  (see  Fig.  1) 
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the  Green’s  function  coefficients  are  expressed  as  follows: 
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Although  (15)  is  expressed  in  integral  form,  its  range  is 
limited  by  the  maximum  and  minimum  values  of  p  so  that  (15) 
can  be  easily  approximated  via  a  polynomial  interpolation  or 
tabulated. 

These  formulas  are  valid  under  the  following  conditions: 

\kD\  >  3fco  \kDp\  >  3.  (16) 

If  these  conditions  are  not  met,  the  lateral  wave  contribution 
diverges.  In  these  isolated  cases,  the  homogeneous  Green’s 
function  for  the  direct  and  reflected  waves  can  be  used  [8]-[10]. 

The  advantage  of  using  this  formulation  rather  than  the 
spectral  representations  [8],  [9]  is  that  the  Green’s  function 


computation  time  for  each  (r,  r')  pair  is  comparable  with  the 
case  of  free  space  since  the  expressions,  although  lengthy,  are 
still  in  closed  form. 
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